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ABSTRACT 


The Coppens-Sanders perturbation solution for the one-dimensional 
non-linear acoustic wave equation with dissipative term deverigune 
the viscous and thermal energy losses encountered in a rigid walled, 
closed tube with large length-to-diameter ratio was extended to include 
sixth order terms. The solution was then investigated to determine 
the region of validity. Computer programs were written to evaluate 
and graph the resulting waveforms. Available experimental results 

eo ee eee 

were compared with the theoretical predictions and good correlation 
was found to exist in the region of low Mach numbers. This agreement 
was found to gradually deteriorate as the Mach number was increased, 


A Fourier synthesis approach is also presented and the leading terms 


of the first ten harmonics are derived. 
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Lagrangian coordinate measured from piston 
(peak) acceleration amplitude of piston 
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Lagrangian coordinate of rigid end of tube 
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forcing term (see Eq. 16) 


G@ = 2H f = (angular) driving frequency 
coe = infinitesimal-amplitude resonance frequency 
AW = GO - ap 
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D for frequency ae 


NX 


operator for body forces (see Eq. 9) 


lems { Sa 
aa 2 C2 Ta = D'Alambertian in Lagrangian coordinates (one 
a c Ot spatial dimension) 


vs = (peak) amplitude of voltage output of accelerometer 
Vu = (peak) amplitude of voltage output of microphone 

Sy = V/A, = accelerometer sensitivity 

Su = Vy/P = microphone sensitivity 
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1. INTRODUCTION 


The thermodynamic equation of state of a fluid for reversible 


acoustic pe van can be written as 
40P) + 5 elf) (1.1) 


or; for an ideal gas, 


jet de oe 
tm (Mp 


(22) 


where J is the ratio of specific heats. The continuity equation for 


one-dimensional wave propagation in Lagrangian coordinates is 


<neen os a Cr 835 


The force equation in the absence of dissipative mechanisms has the 


form ) 
SD ; 
~ ot ir (1.4) 


It is a tede by straightforward combination to obtain the equation 


a § wea Peay . 
-(G)# En Fy /2F oo a 
ae* 


for a fluid obeying Eq. 1. 7 or 
—¢~/ a 
a (4 E oo 
an 


(1.6) 


for an ideal gas. 
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If interest is restricted to acoustic processes for which oA <1, 


Eqs. 1.5 and 1.6 may be approximated by 


= pi DA 2)? Sade 


which (upon differentiation with respect to time) becomes the non- 


dissipative wave equation 


Ther aes a far ve 
a= aeee ba (1.8) 


where 


penotanigel Gy jin Et ae 


ee es 
aa* a Be 


A more general form of Eq. 1.4 is 


3d 2 
fe ae = - a + fk (1.10) 


Beg: is an operator which generates all forces, other than that 


arising from the gradient of the pressure field, which are active in 
the system of interest [4]. If Eq. 1.10 is used instead of Eq. 1.4, 


Eq. 1.7 becomes 


iss yea F shyt Se) ia 


In the Navier-Stokes case 


ig: 
4 = (FUE) fe an 


where 7 and Xp are the shear and bulk viscosity coefficients, respectively. 
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An important contribution to the solution of Eq. 1.11 for the 
case of Eq. 1.12 was made by Fay [5]. He investigated the changes 
in the steepness of the wavefront in a periodic, finite-amplitude 
plane wave of infinite extent in a viscous medium. By means of 
Fourier analysis techniques he was able to show that the nonlinearity 
of the pressure density relationship results in accumulating distortion 
of the waveform with shifts of energy from the lower to the higher 
frequency components of the wave. Because the classical Navier-Stokes 
absorption coefficient is proportional to the square of the frequency, 
the higher frequency components are attenuated more rapidly. Fay's 
solution is valid in a region far from the source where the rates of 
harmonic growth due to nonlinearity and the rates of attenuation due 
to viscosity are tending to balance leading to a waveform of relatively 
stable shape over distances corresponding to many wavelengths. 

Fox and Wallace [1] investigated the same problem as Fay using a 
graphical analysis technique. Their result is essentially the same 
‘as Fay's. Keck and Beyer [6] developed a perturbation analysis technique 
for periodic plane progressive waves of infinite extent in a viscous 
medium with which they generate terms through sixth order. The wave 
equation is the same as that treated by Fay, but the solution is valid 
only near the source. 

Weston [20] has presented a linearized wave equation for the 
propagation of monofrequency sound in tubes, He assumes that the wave- 
fronts are planar except near the walls. The viscous and heat-conduction 
losses in the boundary layer are assumed to be the dominant dissipation 
mechanisms in the fluid. Except in this boundary layer, the wavefronts 


are independent of a radial coordinate. This suggests that the boundary 
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layer loss can be replaced by an equivalent absorptive process active 
throughout the volume of the cavity [4]. More will be said about this 
later. 

Saenger and Hudson [14] present a simple description of periodic 
shocks at resonance using a two part solution. Their solution is 
based on the assumption of the applicability of the linearized acoustic 
equations everywhere except in the region of the shock where the Rankine- 
Hugoniot shock conditions are assumed. They consider the case in which 
the piston oscillates at the fundamental frequency. Compressive viscosity 
is ignored but the effects of shear viscosity and heat conduction in the 
boundary layer are considered. Their solution remains finite only because 
of these considerations. 

Betchov [13] also postulates the existence of the shock wave and 
constructs a solution at resonance based on a continuous and a discontinuous 
part. For an inviscid fluid the amplitude at resonance is found to be 
finite and to be determined by nonlinear effects. The effect of wall 
friction is discussed and it is suggested that this could modify the 
solution significantly. 

Weiss [7] has applied finite-difference techniques to the nonlinear 
inviscid acoustic equations and has shown that repeated reflections at 
the rigid boundaries tend to promote the development of a discontinuity 
in the velocity profile. The absence of any dissipative mechanisms, 
however, eliminates the possibility of any steady-state standing wave 
patterns. 

While most investigations have been confined to traveling waves, 


a few have considered standing waves. 
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Solutions of the nonlinear wave equation have been obtained for 
the case of finite-amplitude standing waves in a closed tube, where 
the tube is driven at one end by a piston source. This work was done 
by Keller [21] using the Lagrangian formulation and assigning an 
unrealistic value of -1 to the adiabatic exponent ¥Y . Due to neglect 
of any dissipative mechanisms, his expression for the particle velocity 
becomes infinite if the piston frequency equals a natural frequency 
of the tube. 

A detailed theoretical analysis based on the nonlinear acoustic 
equations including the effects of compressive viscosity and shear 
viscosity in the boundary layer has been done by Chester [8]. He 
successfully predicted asymmetries in the finite-amplitude standing 
waves. Unfortunately his results are not easily compared to experi- 
mental results. 

Coppens and Sanders [4] developed an extension of the Keck-Beyer 
perturbation approach wherein wall losses as well as bulk losses are 
included. Their solution is applied to the case of finite-amplitudes 
standing waves in rigid walled cavities. The results yield information 
concerning the amplitudes and phases of the Fourier components of the 
waveform, and indicate the importance of the type of absorptive process 
on the resulting waveform. 

The purpose of this thesis is to extend the Coppens-Sanders 
perturbation approach. In part this is accomplished by a direct 
application of their method to obtain sixth-order terms. Information 
concerning hia response of the cylindrical cavity to the sinusoidal 


excitation is extracted from this solution by some algebraic manipulations. 
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A Fourier synthesis approach to the problem is presented and from it 
the leading terms of the first ten harmonics are derived. Information 
similar to that obtained from the perturbation approach is acquired 


and both results are compared to the experimental results obtained 


by Beech [19]. 
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2. THEORETICAL DEVELOPMENT 


In this section both a perturbation approach and a Fourier-synthesis 
approach to the solution of the one-dimensional non-linear, dissipative 
acoustic wave equation are presented. These approaches are formulated 
for a rigid-walled, closed, cylindrical cavity with large length-to- 
diameter ratio. The solutions obtained are for finite-amplitude standing 
waves. Manipulations are also presented which enable each approach to 
yield resonance response information which is convenient for comparison 


with experimental results. 


Perturbation Approach 


Assume a perturbation series Se 2 Sn such that the n-th 
a) 
term is of order n-1l in the Mach number. Substituting this into Eq. 1.8 


and collecting terms of equal order yields the set of equations 


bled D) ter = D, ; Ct) 


where 
ng Ot’ DSy-6 
Pie peg Ow ) jie 
Dp. = 2929 
O n< | 
If Uy= 2s Xnj , where the summation over j is understood to 
d 


encompass only those frequencies included in each Dm, Eq. 2.1 becomes 


2.(O, 4 Dj ) ng = Dy (2.3) 
¢ 
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In order to utilize this formulation, it is first necessary 
to obtain the form of Dj. For the case of finite amplitude sound 
in a duct, Weston [20] has presented a wave equation valid for the 
propagation of monofrequency sound, Weston's equation may be written 


in Eulerian form as 


JY 
Cig dae Geeks = =O ioe 
where 
-/o| |! -| gil 
5+ (3)(2a) ‘hy a4 (1-¥ )y fc (2.5) 
and 


(uniform) cross-sectional area 


perimeter of area S 


(angular) frequency of the propagating waveform 


<7 CSG 
i] 


= kinematic shear viscosity of the fluid 
>. thermometric conductivity of the walls 
irae Es 

This may be rewritten as 

a4 _¢ ou 

ear aso 


When dissipation processes are included, it is usually assumed, explicitly 


ade J “4 
Cet (6 =O aS 


or tacitly, that they are weak. Under this assumption, and for small 
Mach number, the derivatives (Yat and ( Yr yt may be interchanged, 
This also holds for the pair (Sp ya and (Mt)x. Thus Eq. 2.3 must be 
solved with 
(4) ee A ee 
Py «(a oat: ‘Py (2.7) 
where di represents the value of DO for a frequency a; ; 
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The final result is the approximate wave equation 
. 3 
2 Zz 
Fe (ile) SehoBheey ahs - 55 2) any 
: e 
| d da / 
(2.8) 
aot 4> Ja Oa 
Eq. 2.8 was solved [4] for a uniform tube of length L excited 
at end a = O by a piston driven with constant acceleration A, ata 


frequency @) , and terminated at a = L by a rigid cap. In this 


solution, the boundary conditions were given by 


u, (0,0) e (A%Ka) sin (wt am 6; ) 
(2.9) 
Uy Zt) = uy; (9t) = dn; (Lt) =O ,n>l 
where io will be defined later. 


The lowest order solution of Eq. 2.8 is 


U,= Uy, = Ge (uu) (2.10) 


where 


A, of et &,) wf aed Ce if G-a} 


U, ———— ——_——_ C011) 
~ 4 ik “02 L 
er ni <2 
and k = k -i® is the complex propagation constant. 
The solution for wi may be rewritten as 
(2.02) 


+] tinh «(L-a) cca keL-a) | di. wt | 


where 


U, = (2) (ach oa * TNE (2.13) 
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and 


\ 


es ee a Bu bl 
taphk VL 


Eqs. 2.12 and 2.13 will be rewritten in terms of wy, , the 


(2.14) 


resonant frequency which maximizes the pressure amplitude Pu ‘ 
associated with U 


11? at the rigid end of the cavity. 
2 

Oe MM En iss 

‘es Z ma 

wo = wp + Aa 


(25.115) 


The integer m represents the normal mode most strongly excited by the 


. Ge 
input frequency. For oL<< { and aw<«7Uf and ©, are 


= pet 
ae ee 


Cow Ge, 
MT d, Oy 


(2.16) 


0; 


i\s 


Tans’ (~ — }: C25, Lea 
1 Wp 


From this point on, higher-order solutions are obtained by first 
making the approximation 


a ne 
Ui, 


diw K (L-a) Cov wet 


(2.18) 
which is Eq. 2.12 with terms of Wp omitted. The procedure (presented 


in Appendix A), is then merely an iteration process. 


In view of later discussions it is advantageous to point out that 


Eq. 2.18 may be rewritten in terms of the associated acoustic pressure 
by merely making the substitutions 


o > i apes Co 
I a Cet 


(2.19) 
22 


This is equivalent to ignoring terms of order Vp im Bq. L.LO, and 
converts the calculated variable from velocity to pressure, a more 
commonly measured quantity. 

In terms of the pressure, the solution presented by Coppens and 


Sanders is 
tp = ge Kea) ai wt FY G2), (2.20) 
+4 (ZE)*[422.9)-£ (43-2), 
Ui 22eg 42-2] e~ } 
toes 2K) 4 (FP) eh - (Yea 9, 
4(422h]t-} 
2 ¢/, 
too 31-0) 3(F) (23h fH | fe334, 
tf (22354), + (223,45), 
+3 (2,23), | 4 a ;, 
paca Kea QP [ia 54huty C24, ]r~-} 
toa 5K U-a) f(A] (23.45), 
+g 224,55 +5 (82 55)s [4 ; 


hae 
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where 


(4, 4)s Mes a din (gut + G+) (2.21) 
- ree 
Po 


oe =bf. 7 
G6, a Cot ee ON oe Zauy, /; Wey 2 8 aay 


On /a. (2.24) 


and M is the peak Mach number of the first order solution. 
Using the same scheme, the sixth order terms were derived. 


The solution for Pe which is to be added to Eq. 2.20 is 


ie = Con 2KLa aH) l22ssth i255 


24 


Ch - 
oe Ah Y Keo) 3 (BY /e Z44¢,5)y 
47 ( 33445 )4 ta 545). 
+(Z5344)9 + £(223,44)¢ 
ide 7 bdland (heck Jy 


42 (42259 )4 +2 (22234) y 


/ 
+5 (1222¢), 14 (13e3¢h fe] 


or Kika 147 fee 3.4 5.) 


17 22456), +3 (22356), 


, 


PEPE eS ae 3(2224b) 


Pee vz 2 36) Jens 
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The entire perturbation solution herein developed, which is the 
sum of Eqs. 2.20 and 2.25, will be discussed in some detail after 


the Fourier synthesis approach has been presented. 


Fourier Synthesis 


Beginning with the approximate wave equation Eq. 2.8 and assuming 
the total solution is composed of the known classical solution & 


satisfying 


Ty e+4 a Je zt 7 ace 


a 3c ip d fc 


=O (2.26) 


and some other portion ¢° defined by 


iF SMEARS (2.27) 


then 


(2.28) 
- y+ a: 
= 2 sees ere 


Define the phase velocity i to be decd es 5) and Eq. 2.28 may 
be written 


re ae ta ES 
ia 2 ee (SAE 


yi rv) os : ae 
= 9 Jessie a Vail 


(2.29) 


Assume J may be written as a Fourier series 


= “1 Cor nKL- -2) cot (nwt + Pn ) (2.30) 
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and 


& = Xe Cor K-a) Cot wt (2.31) 
OW 
and let 
go. i 
=) =)! Arla) coa(nust + fn ) (2.32) 
hrs 


—T 
where An and |» will be obtained later. Combination of Eqs. 2.29 


through 2.32 yields 
neon 20 fas 2 ay. 0 
es aes bin ™K(L-4) or frat + Br ih 
—dnXn din AKL-&) dia rwt + Zn’) (2.33) 


a i zp An @) cor (nwt +/m ) 


We impose the conditions that 


Ar,(@) = Ba cow nk (L-2) (2.34) 


and 


An'(a) < HK din 2Kb-2) 


so that Eq. 2.33 becomes 


ULC) lon] eabate 9-5, salt) 8 
= 7 Dy Cor (nut + By 


Introduction of the phase angles Gy of Bq..\2.17 allows Eq. 2.35 to be 


(2.35) 


written 


ae pin (reat tg -On) = if ee Bm cor frwt + ) (2.36) 
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: aa 
In order to simplify the manipulations to obtain eh and /y 


define 


X,corK (L-2) cout t $i) =X Cos KL-2) Coslut + 9 


i Xe Cot K(L-a) Cor wt 


Xn Cor nk(L-«) Corfiut t Px ) 


= Xn Cor nk(La) Cor (ruts bn”) ea 


Thus 


$v 
ae Y Xn cor MK La) Con (rt + br ) 


n= 
and 
a Se 
Le xy a er. 4; ) 
de a Xa Xoneg cot (nwt t $y. -j +¢; ) 
ie 
Combining Eqs. 2.40 and 2.36 gives 
- £ me Xie Li.lnut+ ox - 6») 


Pa zi Con (nwt t+ Dany z -47) 


r= 


ae Xj Xn-4 Cos fut + by-j +; ) 
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(2% 


(2 


(2. 


O71) 


-40) 


It is convenient to shift to complex notation at this time. Let 


X, Cor 4, & x, Cot f° + Xe 


(2.42) 
X, bin Ji i J ough er 
so that 
hee 
pie + Xe (2.43) 
and 
= Z ° 
my, gee anu y 2 vn | aay 
Then Eq. 2.41 becomes 
4. yo i (tr-en-%) 
St! Cor On 
(2745) 


ei: i(b-ntj ok 
pthepierton b.D iphag ltrs Hy 
j=l 
If the equivalent pressure equation is desired, make the 


substitution 


a = Yo? (2.46) 


where M is the amplitude of the Mach number for the particular harmonic. 


Define 


OP a 


so that Eq. 2.52 becomes 
° - 
My L Mn On 2) 


Dy 
Se a — Py 
ae ES a (fnej Fj) (2.48) 
Gag age eo Ong 147) 
ep iis aes ee ea 
iv 


Examination of Eq. 2.48 and Eqs. 2.20 and 2.25 shows that the two 


methods are equivalent since expansion of the term 
n-/ 
Z Zo TG igen gr 
yields the iaqeiar aaa of each harmonic derived in the perturbation 
scheme, while the other portion of Eq. 2.48 contains the correction 
terms. 

Recall that in both the perturbation approach and the Fourier 
synthesis approach the leading term of the fundamental was calculated 
to be ue K(L-a) Grngue . This is of order one in the Mach number. 
It is also a known experimental fact that the fundamental remains 
essentially unchanged with respect to Mach number. This suggests that 
the phase angles associated with the corrective terms might cause the 
corrective terms to sum to approximately zero, It then becomes 
conceivable that a solution composed of only leading terms could 
yield significant results, In order to investigate this possibility, 


Eq. 2.48 was used to obtain the leading terms of the first ten harmonics. 


The result is 


% = D Cov Ka) | 2 wt} (2.49) 


+ 3 cor 2K(ka)t (LP ) (2), 
1 MB si (23 
ty vow 4Xtlea) (S [23% 
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% 7 (24 +)u]} 
tf Cor 5K (La) f 4 Sh)'[i24 45/5 
if (224shofeassh]} 
ef tor be (La) (He s45b) 
pee: +2 2350) 
+Z(Z2446) +S (22246), 
Re (2252 od 
$f cco Dea) (SF )'fez 4567), 
tg 224567) + HELLY, 
+3 (223467), + 3 (22846), 
a a, tz (22 3459) 


ad CEIYS1), +4 (222357), 
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#3 (223347), +B Ges 
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4 fect sKLall(F) fi2345 624), 
#7 (2245579) +3 (2235678), 
+2 (22246 7)q + 5 (222967 Ep 
+ £2245 5958p + 5 (22 29578)s 

3 (22299 + J (2235728) 


A ae 


45 @Z53¢ 


BE) esl 

Ye +f (2225478) 
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/ = 
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) Saal 


a ~ / 
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2 (2254561 tp + J (23245679), 

tg (22235,619), + £ (2224679) 
+i (22224 SD )y # S (22223699), 
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where (a,b)g has been given in Eq. 2.21. 
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Cavity Resonance Response 


Investigation of the behavior near resonance of each approach 
leads to the desirability of obtaining the response of the cavity 
for each harmonic as the fundamental is swept through resonance. 
itis obvious from the” form’ of Eq8! 2.2/7, 2.32; and 2.55, that) both 
approaches yield results which may be written as sums of terms 
having like factors of K(L-a) andwWt. It is also evident that each 
of these harmonics is heavily dependent upon M, the Mach number, and 
Aw. 

The algebraic process necessary to bring each of the two solutions 
to a form which conveniently yields the desired information will be 
derived in a general form and then applied to the perturbation 
solution, It is then easily extended to the Fourier synthesis solution. 


We begin with the relationship 


~ 
D4 Lie Got +97) a AL Gut+d). (2.50) 


At t = 0 Eq. 2.50 becomes 
~ QAY 
a : - bi ste e : 
> @é bin $e = Adin $. 
ioe 
The derivative of Eq. 2.50 evaluated at t = 0 yields 
NY 
aoe ¢, a 
La Le Ct G = feor $ (2.52) 
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which when combined with Eq. 2.51, yields 


A= (52a: hi tc) + fas i Cot ¢)* (2.53) 
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The sum of Eqs. 2.20 and 2.25 may be written as 
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ez A coe g tee) bi. Gott 9. 
for each harmonic, where j indicates the harmonic number and n 
indicates the number of terms in that harmonic. It is evident 
that the amplitude of each harmonic may be obtained by applying 
BW. 2053.5 

In applying Eq. 2.53 to the perturbation approach the magnitude 
of each harmonic was normalized to approximately one by the maximum 
magnitude of a factor chosen from the leading term of that harmonic. 
This normalization was done to yield a set of curves whose maxima are 
of order unity. For the fundamental this normalization factor is one. 
The method of choosing the normalization factor for higher harmonics 
will be demonstrated for the second harmonic. The leading term of 


the second harmonic is 
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and the leading term of the second harmonic may be written 
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Since 
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Eq. 2.59 may be written 


Aefh GA Ce G2 Lnfest #82) (2.61) 


2 Tl p Co cy dz 
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and the normalization factor for the second harmonic is then chosen as 


eeien. orse 
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It should be noted that these normalization factors are not the maximum 


(2.63) 


amplitudes of the various harmonics, neither are they the maximum 
amplitudes of the leading terms of each harmonic; rather, they are 
factors chosen to yield maxima of order unity. 

These factors are given the symbol MAXA | where n refers to the n-th 


harmonic so that 


fof th 
TAXA , ae Lm I Wy Ce Br ao. (2.64) 


By combining Eqs. 2.53 and 2.64 it is possible to define a. as the 
normalized response of the n-th harmonic. It should be noted that A 
defined in Eq. 2.53 is a function of the Mach number and AW, 


Thus 


Q, (M Aw) = len (2.65) 


where n is the harmonic number. 
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3. APPLICATIONS 


Associated with each of the solutions formulated in the previous 
section are three computer programs. These programs are concerned 
with the pressure waveform, the Q-curves, and the phase information. 

The perturbation solution has the programs FINAMPI, QCURVES, and PHAMP 
connected with it while FINAMPIZ, QCURC and PHAMPC are related to the 
Fourier synthesis. 

The first attempted utilization of the perturbation solution 
consisted of writing a computer program involving all terms through 
sixth order. This program, entitled FINAMPI, was designed to compute 
and graph both velocity and pressure waveforms. FINAMPI was first used 
to graph the waveforms contained in Ref. [4] and has been used since 
in determining for what Mach numbers the theory begins to deviate 
significantly from the experimental results. The graphical output 
from the computer is contained in Figs. C.1 through C.4. 

The program designed to compute and graph the cavity resonance 
response for the perturbation peiat ven is entitled QCURVES. This program 
is a direct application of Eq. 2.65 and has been used to obtain the 
curves contained in Figs. C.5 through C.22. 

PHAMP, the third program associated with the perturbation solution, 
was designed to compute the phase angles of each harmonic as the 
fundamental was swept through resonance. It is a direct application 
of Eq. 2.54 and resulted in the data contained in Table 4.2. 

In an attempt both to determine the importance of the corrective 
factors to each harmonic and to test the validity of the Fourier synthesis, 
computer programs involving only the leading terms of the various harmonics 


were written, Unfortunately a program involving all ten harmonics could 


Al 


not be brought to fruition because of time limitations. A second 
program, FINAMPIZ, which uses the leading terms of the first six 
harmonics has been used to graph the pressure waveforms also contained 
in Figs. C.1 through C.4, 

Eqs. 2.65 and 2.54 were aie applied to the Fourier synthesis 
yielding the programs QCURC and PHAMPC. QCURC computes and graphs 
the Q-curves for each of the harmonics while PHAMPC computes and graphs 
the phase dependence of the individual harmonics on AW. 

It should be noted that of the two programs concerned with phase 
dependence, only PHAMPC graphs its output. This was done because the 
Q-curve programs indicated the Fourier synthesis approach produced 
curves which were in much better agreement with the experimental results 
than did the perturbation approach. The graphical ae PHAMPC is 
contained in Fig. C.23. 

In this section only a superficial description of each program 
was given. A detailed description of all programs is contained in 
Appendix B, and a more complete discussion of results follows in 


Section 4. 
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4. RESULTS AND COMPARISONS 


In this section the theoretical predictions obtained from the 
applications described in Section 3 will be compared to the experimental 
results obtained by Beech [19]. 

Before proceeding to a discussion of these results it is necessary 
to attempt some qualitative evaluation of the apparent properties of 
the perturbation solution. The amplitudes of the (n + 1)th order terms 


tend to behave as [4]. 
a 
n ~\ 
a) a Mp “Al oe 
” ieee E ge ae (4.1) 
Z An / r= Z 
iré / 
It is clear that this expression can diverge for increasing n. This 
divergence, which occurs because of the dependence Ins Vow, would 
not occur if dy, >d, since the left hand side of Eq. 4.1 would then 
n 
be bounded by (AS) and all series would converge for Mp< Be 
? . Ba: : M 2 
The Mach number which satisfies the expression Zh / and would 
therefore appear to be related to shock formation is M = 0.02. This is 
within the same order of magnitude as the experimentally measured Mach 
number at shock formation, M = 0.01. 
The divergence may be remedied if it is recalled that the bulk 


absorptive processes have been neglected up to this point. If the 


viscous terms are retained and assumed to be additive to the wall losses 


then 
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The analysis and results are the same as before with the H's and 9's 


redefined as 


dnt uve 
fh epee (4.4) 


(4.5) 
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With these modifications the amplitudes of the (n + 1)th order terms 


now behave, for sufficiently large harmonic, as 


ht or 


which tends to zero with increasing n. 

The modifications introduced by Eqs. 4.2 and 4.3 are important 
in the perturbation expansion for those harmonics whose indices 
satisfy the inequality ne Ye . For the problem dealt with 
in this research the theory developed without the bulk absorptive 
terms is adequate as long as attention is restricted to frequencies 
below about the thousandth overtone. 

In comparing the theoretical and the experimental results it is 
convenient to divide the discussion into three areas: the low Mach 
number region, the intermediate Mach number region and the high Mach 
number region. 

In the low Mach number region, below M~0.006, it is obvious 
from Figs. C.1 and C.2 that theory and experiment are in very close 
agreement. In these figures of pressure waveforms the large drawing 
contains the waveforms predicted from the computer programs and the 
smaller inset is a drawing from an oscilloscope photograph of the 
waveform observed experimentally. The dashed line in both drawings 


is a sine wave which is included for reference purposes. One of the 
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main reasons for the good agreement in this region is due to the 
fact that the ("P/25, ) terms are as yet relatively small. 

It should be noted here that in all the figures of pressure 
waveforms, the major difference between the FINAMPI and FINAMPIZ 
predictions is in the phase dependence. Recall that the development 
of the reasoning behind the FINAMPIZ program rested on the possibility 
that corrective terms to the generated harmonics could be ignored. 
These curves lend support to that suggestion in the low Mach number 
region. 

The low Mach number region being discussed also includes the 
first two sets of Q-curves, Figs. C.5 through C.16. It is evident 
in both sets that the agreement of the experimental results with the 
predictions of the Fourier synthesis approach is excellent in the 
lower harmonics and deviates only slightly in the higher harmonics. 
When the agreement does deviate, the experimental results in general 
fall between the two theoretical predictions. This tends to indicate 
that the two solutions have the exact solution bracketed. It is 
evident in the QCURVES predictions for M = 0.005 that the third 
harmonic curve is no longer a smooth bell-shaped curve. The reason 
for this is the growth in the third harmonic of the correction terms 
arising from the fifth-order solution, For instance, the ratio of 
the approximate magnitude of the first corrective term to the leading 
term in the third harmonic is about 0.5 at M = 0,005. This is no 
longer insignificant so that the predictions based on the perturbation 


approach begin to show unrealistic fluctuations. 
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Theoretically predicted amplitude information is obtained from 
the two programs, PHAMP and PHAMPC, The amplitude results from these 
programs together with the experimental results are contained in 
Table 4,1, It should be noted that Beech's uncertainty of + 0.1 Hz 
in frequency results in an uncertainty of + 0.628 radians in Figs. C.5 
through C.23. A comparison of the uncertainty in amplitude (based 
on the uncertainty in 4W) with the observed discrepancy between 


measurement and theory for M = 0.004 and4W= 0 yields 


n uncertainty in discrepancy 
harmonic content between PHAMPC 
and experiment 


2 OFS 0.4 

J O53 Ot? 
4 OTS (Oa be 
5 0.07 0.08 
6 0...03 0°. 07 


It is evident that the uncertainties in harmonic content are consistent 
with the discrepancies in harmonic content noted between PHAMPC and 
experiment. 

As was stated in Section 3, PHAMPC also predicts the phase 
dependence of the various harmonics. The results of this program 
are contained in Table 4.2 and Fig. C.23. One thing should be noted. 
The sets of phase angles predicted by PHAMPC are identical for all 
Mach numbers, which is not surprising if the reader recalls that the 
normalization process destroyed the Mach-number dependence of the 


leading terms of each harmonic. 
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One other program exists for the investigation of phase and 
amplitude characteristics. The computer program FOUANAL is being 
introduced here since it was not derived from either theoretical 
approach. FOUANAL takes data from an oscilloscope photograph and 
computes the Fourier coefficients of the various terms of a Fourier 
coefficients of the various terms of a Fourier expansion of the 
waveform. The results from this program are also contained in 
Tables 4,1 and 4.2. 

It is clear that the results from FOUANAL are fairly accurate 
in the lower harmonics (first, second and third) but tend to become 
much worse in the higher harmonics (fourth, fifth and sixth). Judging 
from the amplitude results it does not seem likely that FOUANAL would 
yield accurate results in the phase angles of these higher harmonics. 
Examination of Table 4.2 verifies this supposition. These erratic 
results in the FOUANAL calculations may possibly be explained in the 
following manner. The input data for the program come from extracting 
sixty-four points from a standard oscilloscope photograph. This was 
done by hand and was subject to many errors. It is possible that the 
data were not given to sufficient accuracy to enable the computer to 
predict accurately past the third harmonic. The most feasible suggestion 
for improvement in this respect is to use a digital counting procedure 
to obtain the data. It should be noted that the data in Fig. C.23 
are the FOUANAL results for M = 0.004. Clearly these points are 
extremely random and the uncertainty in the experimental MW cannot 
possibly account for the entire discrepancy. 

Shifting now to the medium Mach number region (0.006< M< 0.009) 
we see from Fig. C.3, for M = 0.007, that although the general agreement 


is still good, deviations from the previous good agreement are becoming 
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evident in the computer predictions. Unfortunately experimental 
data for Q-curve and phase angle comparisons were not obtained 

but it is not unreasonable to assume that these comparisons would 
show essentially the same result as stated for the pressure waveform. 

Moving into the region of high Mach numbers, M = 0.009, we find 
that the pressure waveforms predicted by each program are generally 
invalid. This is nearing the experimentally observed Mach number for 
shock formation (M = 0.01), and since the theory is predicted on small 
Mach number, pre-shock conditions it is not surprising to see these 
deviations, 

Essentially the same results are observed in the Q-curve comparisons, 
Figs. C.17 through C.22. The reason for the experimental amplitudes 
falling off so rapidly is at present unknown. Tables 4,1 and 4,2 
again reflect the serious discrepancies in the high Mach number region. 

It is now possible to make some general comments concerning these 
two approaches. Both the perturbation approach and che Fourier synthesis 
approach have adequately described the effects of finite-amplitude standing 
waves in rigid-walled cavities provided the Mach number is relatively low. 
This precludes that area in the pre-shock region where significant 
distortion exists. It has been demonstrated that in the region where 
the approaches are fairly accurate, Q-curve information can be easily 
obtained. This information is convenient for comparisons with experi- 
mental results. It is also possible to obtain some information on the 


phase angles of the harmonics. 
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It seems plausible now to suggest the following: Since the 
Fourier synthesis approach gave the most accurate information, the 
computer program involving all ten harmonics should be reviewed in 
the hope that a usable program will result. This program could be 
used to investigate the effects of the higher order harmonics on 
the existing solution and possibly give an indication of whether 
extending the present developments past six harmonics is beneficial. 

Finally, after some exposure to basic communications theory, it 
appears that a good grounding in this subject is mandatory if the 
present theoretical line is to be pursued further. It would not be 
surprising to find that a more exact solution may be obtained as a 


result of the techniques used in communications theory. 
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APPENDIX A 


Coppens-Sanders Iteration Scheme 


From Section 2, the approximate wave equation for finite-amplitude 
standing waves in a ite walled tube is 
a Ca + G es Bee : ) dng 
uf antot G da 


7 (A.1) 
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For convenience the rar hand side of Eq. A.1 will be written 
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i are complex, and the three dots in the 


exponential term represent phase terms to be defined. Work is done 


where the ay and 


in one n at a time so that the n subscript may be dropped. 


Expansion of the left hand side of Eq. A.1 yields 
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where (-G) was previously defined to be Gr.) 


Thus Eq. A.1 becomes 
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for a given order n. If Eq. A.5 is generalized in time it becomes 
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Substitution of Eq. A.3 into Eq. A.6 yields 

Pages: Jie ee ip & 2 
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Now each Y will be of the form 
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Substitution yields 
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Additional substitution of Eq. A.9 and Eq. A.3 into Eq. A.7 gives 
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Proper manipulation of the bracketed term in Eq. A.10 will yield the 


(A.9) 


phase relationships mentioned previously but left undefined, 


First recall that 
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(A.14) 
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Substitution into Eq. A.10 yields 
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It is now possible to write down the now U as 
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Define the following quantities; 
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Tow iy = i (A.18) 
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Pictorial representation of the quantities 
defined in Eqs. A.17 through A.19. 


Figure A.1 
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With these definitions, Eq. A.16 becomes 
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and 
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yielding 
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The preceding was a general outline of the procedure to be 
followed in the Coppens-Sanders iteration scheme. The remainder of 
this section contains a derivation of the terms through third order. 


Assume the input 


Uy = ff “he Mea) Coa wl (A.23) 


which was shown in Section 2 to be the calculated first-order solution. 


Then we have 
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The next step is to obtain x2° 


/ 
Sr OK n-¢ ' 
‘2 ae ‘ d o $; ee Os: os, 


ag da Ja da (A.27) 
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Substitution of Eq? A, 26 IntevEq, A.2/ yields 


Lz = /7 ee eee SO om 


(A.28) 
where M = An /Co 
Expansion of the trigonometric term gives 
Cor ?> Kk-2) tu.*wTt (A.29) 


= CU feu fates} tiffer(v)-1] 


from which the constants + 1 and - 1 may be discarded since application 


2 
of Yadt will remove them. Then 


Z =). soe, Cor MC-«) Coe LW (A. 


which implies 


A22 = - ¢9* (A. 


and from Eq. A.20 
ce, ea ee C62 


Zo ae 


Recall that Mag @ 4n7 , and hence, 


> UR ae Ls Kh -a) Cor (Zu +L) (A. 


which is the second order solution. 


The general equation for Iq is 


<i == Sen m Mea) tu (mutt). (a. 
Substitution of Eq. A.34 into Eq. A.21 results in trigonometric 
such as 
Cer. mk La) di.fnt by) Coe Nk a) tir lat +h; ) (A 
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(A. 


30) 


31) 


32) 


33) 


34) 


products 


ae =) 


For simplicity let 
X = K(h-a) 
rs wt (A. 36) 
Then Eq. A.35 may be written 
Coa mX Li MT + dy) Cor HX te. (C744) (A.37) 
Expansion of Eq. A.37 gives 


Cf foo (m4 )X Cor (04) Tt byt Ge} 


(A.38) 


— Cor (nt) X Got rT + he by f 


where terms of the form sin mk(L-a) cos¥VC have been ignored 

in Eq. A.38 since substitution of these terms into the perturbed wave 
equation result in amplitudes of order Ye which have been neglected 
throughout. Eq. A.38 must always be written as (m - 1) where m > ‘ 
since a positive frequency has always been assumed, 

Before proceding to the derivation of the third harmonic terms, 
a brief discussion might be of some use in acquainting the reader with 
higher order terms. 

In the third harmonic derivation, the user of this iteration 
scheme is first exposed to the development of more than one term in 
the solution for a particular order. In other words, the iteration 
yields both the leading term to the third harmonic plus the first 
correction term to the present solution. At the same time, the use 
of Eq. A.38 is demonstrated for the first time. It should also be 
pointed out that the methods used in the third harmonic iteration are 
identical to those used in the higher order iterations, regardless of 


the number of terms encountered. 
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Following the step by step procedure then, 


2 
Par ts A ISz-¢ Di JS OR oF = Os gh 
at ety ‘Tr. 7 ie ae Bees ee 


Grr, 
Remember , 


--& ten Kl-a) bin, wt 
&/) 


ai = ee cor Ka) abi. WE 
i= os (SF) 7 Pdi. 2-4) de (Zat + #) 


oe BE eer, cor 2K C-a) yl Ct G2) 
ae UNG AE ae 


so that 


Wi 72(-# - Ces K(L-a) divatl Ld 


noe 2K) bite (et éz)) 


or using Eq. A.38, 


‘ai af & a5 [fon PR(L-4) Caf Fut + G2) 


— Cor K(<-2) cme (utr 62) | 
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-40) 


-41) 


42) 
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Using Eq. A.8 implies 


£ OP WF 
fae ke Las 


heile wen (A.44) 


/ 
eo 2 
aa a (A.45) 


and from. Eq. A.10, 


Uae = y L th We 2 ae fees AS (A.46) 


Co CG. 


= Ub ae OY 


and 


Co 


go hOt82) 
2 Tl. oe dite © 


Now Uz = Uz; tX3} , and when the substitutions are made, the 


Us, =~ z Us Ga 2 Un oy fp & e 


(A.47) 


final form of the third harmonic iteration becomes 


Upc 3 Lp wr hh IK Ea) 0orCul 4B +63) 


(A.48) 


~ LU (Yyoh ds Hien) colt 4s¥6s). 
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A combination of the first three orders to obtain the total 


Solutuon thus fan results an 


“= Qh iw KL-a.) Cn we 
Af, ML) di 2K J y, (A.49) 
tn (G A = coe (Lu 462) 


7 Ot) ey ti. Skb-4) co Cat 4B2 463) 
- £U (4) Ty i. K¢-a) Cor(al + G+ 62). 


Notice that the second term of the third order solution is of order 

one in j, (with reference to Eq. A.3), and thus it is termed the 

first correction term to the fundamental or first order solution. 
Rewriting Eq. A.49 in order to achieve a form consistent with 


this last observation and the form of Eq. 2.20 results in 


“9 " * dw Lika) Lee at - 


Si (4) coful 46,7 62) 
Fi .f 


+ di. 2a) bd (BE) fom (te &2) 


pe) 


4 di. MESS (GE) Ty, Cor (Sat + G2 + @) 
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(A.50) 
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APPENDIX B 


Computer Programs 


An IBM-360 digital computer was utilized in performing the 
calculations, and obtaining the wave forms predicted by the theory. 
Computer programs were written for each of the basic formulations 
as well as for each of the variations investigated. Each of these 
programs, together with a short description, is presented in this 
section. Included with the description is a list of the symbology 
peculiar to that particular program. 

B-1: Program FINAMPI 

This program has written to perform the calculations and plot 
the graphs for the perturbation solution. It was originally written 
in FORTRAN 60 language and subsequently converted when the IBM-360 
became available. 

The program was designed so that the Mach number, resonance 
frequency, length of the tube, position in the tube, 7, Aa), and a 
are the input parameters. The first half of the program is devoted 
to the velocity profiles and the second half to the pressure waveforms. 
The functions defined at the end of the main program are part of the 
conversion from FORTRAN 60 to FORTRAN IV. A list of the symbology 
pertinent to this program, together with the corresponding parameter 


in the original formulation follows. 
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Table B.1l 


D(N) - Sm 
S = Mach number 
WD = Aw 
WR = Wp 
: ee 
F = L 
A =e 
U = Velocity 
P = Pressure 
PL = © 
WT _ wt 
as) 
OL)» i= aie s /2Hg 
THET(N) = Gy 
Vxyz = the expansion of (X,y), for the velocity equation 
Yxyz = the expansion of (X,y), for the pressure equation 
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B-2: Program FINAMPIZ 


This program was designed to duplicate the output of FINAMPI 
using the Fourier synthesis formulation through sixth order terms. 
Since the final results of the perturbation approach are in the same 
form as those of the Fourier synthesis approach the program FINAMPI 
was converted to FINAMPIZ by setting the corrective terms to zero. 


Because FINAMPIZ is essentially FINAMPI, Table B.1 is applicable here. 
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B-3: Program QCURVES 


This program was designed to compute and graph the normalized 


Q-curve for each of the six harmonics using the perturbation solution. 


In this program, the input parameters consist of the resonant 


frequency, the tube length, JS, P, and the piston acceleration. The 


symbology pertinent 


AO 

F 
Delt(N) 
H(N) 
THET(N) 


SCRIPN(M) 


to this program is given in the following table. 
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length of tube 
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hypothenuse of triangle described in Fig. A.1l 
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B-4: Program QCURC 
This program is the Q-curve program based upon the Fourier synthesis 
approach using only the first six harmonics. The curves are normalized 


in exactly the same manner as in QCURVES and no new quantities have been 


defined. 
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B-5: Program PHAMP 


This program was designed to compute the amplitudes and associated 
phase angles for each harmonic. Essentially this program is the QCURVES 
program with only minor modifications. One new quantity has been defined 
as PHASE(I,J). This is the phase angle associated with the i-th harmonic 


and the j-th Aas, In all other respects Table B.2 is applicable. 
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Liat 16,50 sF 13 210s Fl os 1 


COMPUTES THE AMPLITUDE AND PHASE OF 


PORTION OF THE PROGRAM 
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B-6: Program PHAMPC 


This program is identical to PHAMP except that it was written for 
the Fourier synthesis approach and graphs the phase angles associated 
with each harmonic. In order to use the graphing subroutine the 
subscripts on the phase angle had to be reversed. PHASE(I,J) is now 
defined as the phase angle associated with the i-th Aw and the j-th 


harmonic. In all other respects Table B.2 is applicable. 
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B-7: Program FOUANAL 


This program is an adaptation of work done by Roy M. Johnson, 
Assistant Professor, Naval Postgraduate School. Basically the program 
takes a waveform obtained directly from the cavity and performs a 
Fourier analysis on it. The output of interest is a set of Fourier 
coefficients and their associated phase angles. The program is 
written in such a manner that a table of definitions is not appropriate. 


Each variable is defined within the program. 
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Q-curves for the fundamental, M = 0.004, ~ 
QCURC prediction,---QCURVES prediction 
Experimental results are indicated by @, 


Figure C.5 
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Q-curves- for the second harmonic, M = 0.004 
——QCURC prediction,----QCURVES prediction 
Experimental results are indicated by @ 


Figure C.6 
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Q-curves for the third harmonic, M = 0.004, : 


QCURC prediction,----QCURVES prediction 
Experimental results are indicated by @. 


Figure C.7 
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Q-curves for the fourth harmonic, M = 0.004, 


QCURC prediction,----QCURVES prediction 
Experimental results are indicated by @- 


Figure C.8 
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Q-curves for the fifth harmonic, M = 0.004, 


QCURC and QCURVES prediction 
Experimental results are indicated by @. 


Figure es 
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Q-curves -for the Sixth harmonic, M = 0.004 


QCURC and QCURVES prediction 
Experimental results are indicated by@, 


Figure C.10 
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Q-curves for the fundamental, M = 0.005, . 
QCURC prediction,----QCURVES prediction 
Experimental results are indicated by@. 


Figure C.11. 
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Q-curves for the second harmonic; M = 0.005, 
QCURC prediction,----QCURVES prediction 
Experimental results are indicated by @ — 


Figure C.12. 


134 


& 
@ 
Ro} 
® 
= 
As ‘s 
‘ 
% 
e ‘ 
\ 
A 
\ 
J ki 
iy ‘=. 
7. 
/ . 
x 
/ ‘ 4 \ 
/ N 
/ ss 
/ \ 


-2 QO 2 


Q-curves for the third harmonic, M = 0.005,. 
QCURC prediction,---QCURVES prediction 
Experimental results are indicated by @ 
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Figure C.13 
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Q-curves for the fourth harmonic, M = 0.005, 


QCURC prediction,---QCURVES prediction 
Experimental results are indicated by @ 
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Figure C.14 
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Q-curve for the fifth harmonic, M = 0.005, | 
QCURC and QCURVES prediction 
Experimental results are indicated by @.. 
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Figure C.15 
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Q-curve for the sixth harmonic, M = 0.005, , 


QCURC and QCURVES prediction 
Experimental results are indicated by @e. 


Figure C.16 
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Q-curves for the fundamental, M = 0.009. 
: QCURC prediction,---QCURVES prediction 
Experimental results are indicated by@. 
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Figure C.17 
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Q-curves for the second harmonic, M = 0.009, 
QCURC prediction,----QCURVES ‘prediction 
Experimental results are indicated by @. 
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Figure C.18 
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Q-curves for the third harmonic, M = 0.009, 


QCURC prediction,----QCURVES prediction, 
Experimental results are indicated by@, 


Figure C.19 
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Q-curves for the fourth harmonic, M * 0.009, 
QCURC prediction,---QCURVES prediction 
Experimental results are indicated by @, 


Figure C.20 
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Q-curve for the fifth harmonic, M = 0.009, 
QCURC and QCURVES prediction 
Experimental results are indicated by @. 


Figure C.21 
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Q-curve for the sixth harmonic, M = 0.009, 
QCURC and QCURVES prediction 
Experimental results are indicated by @, 


Figure C.22 
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Phase dependence of the various harmonics on A@® 
as predicted by PHAMPC. n = harmonic number, @ 
= data from FOUANAL for n'-th harmonic 
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